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We present a general expression for the association rate for partially diffusion-controlled reactions 
between spherical molecules with an asymmetric reactive patch on each surface. Reaction can occur 
only if the two patches are in contact and properly aligned to within specified angular tolerances. 
This extends and generalizes previous approaches that considered only axially symmetric patches; 
the earlier solutions are shown to be limiting cases of our general expression. Previous numeri- 
cal results on the rate of protein-protein association with high steric specificity are in very good 
agreement with the value computed from our analytic expression. Using the new expression, we 
investigate the influence of orientational constraints on the rate constant. We find that for angular 
constraints of ~ 5°-15°, a typical range for example in the case of protein-protein interactions, the 
reaction rate is about 2 to 3 orders of magnitude higher than expected from a simple geometric 
model. 
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I. INTRODUCTION 

The association of two macromolecules, in particular 
the formation of protein-protein complexes, is an ubiq- 
uitous process in biology. In the simplest case of the 
associating species being modeled as uniformly reactive 
spheres, the diffusion-controlled association rate is given 
by the classic Smoluchowski result lj, fcoc — 47r£)i2, 
where D is the relative translational diffusion constant 
and R denotes the sum of the radii of the molecules. 
Typically, however, successful complex formation hinges 
on the proper relative orientation of the reactants, which 
can be represented by molecules carrying reactive surface 
patches that have to come into contact with high steric 
specificity for the reaction to occur. 

The simple approach of multiplying the Smoluchowski 
rate constant for uniformly reactive molecules by the 
probability that in a random encounter the two molecules 
are properly oriented ( "geometric rate" ) yields rate con- 
stants that are commonly several orders of magnitude 
lower than the observed values. Some authors attributed 
this puzzling behavior to the presence of long-range at- 
tractive interactions between the molecules that not only 
generally speed up the rate of encounter of the molecules 
but also help "guide" the molecules into configurations 
close to the proper mutual orientation. 
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In addition to this approach, various attempts have 
been made to quantitatively elucidate the influence of 
orientational constraints and rotational diffusion on the 
association rate constant. Among the earliest studies, 
Sole and Stockmayer derived a formal solution 0] of the 
association rate constant of spherical molecules with ax- 
ially symmetric distributions of reactivity and presented 
numerical results for the simplified case of one of 
the molecules being uniformly reactive. Schmitz and 
Schurr Q investigated both analytically and numeri- 
cally the problem of the reaction between mobile ori- 
cntable spheres, carrying single, axially symmetric reac- 
tive patches on their surface, with localized hemispheri- 
cal sites on a plane. Shoup et al. [j| introduced a gen- 
erally applicably approximative treatment that allowed 
simplification of the complex formal solutions of Sole and 
Stockmayer 0] and Schmitz and Schurr Q to closed an- 
alytical expressions; this approximation was also used 
by Zhou [fj in deriving an expression for the association 
rate when each molecule bears an axially symmetric re- 
active patch. All these approaches showed that, because 
of relative angular reorientations caused by translational 
and rotational diffusion, the reduction in association rate 
brought about by orientational constraints is significantly 
less than suggested by the reduction in the probability for 
a properly oriented encounter. 

The previous analytical treatments, however, impose 
only (at most) axially symmetric orientational con- 
straints, whereas no analytical treatment has been pre- 
sented thus far for the general case of asymmetric reactive 
patches (as in the important case of sterically highly spe- 
cific protein-protein interactions), where the precise rela- 
tive orientation of the binding partners has to be specified 
and appropriately constrained. 
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The only numerical estimates for the association rate 
constant for this general case stem from Brownian 
Dynamics simulations, as for example performed by 
Northrup and Erickson who consider diffusional as- 
sociation of spherical molecules, each bearing a reactive 
patch composed of four contact points in a square ar- 
rangement on a plane tangential to the surface of the 
molecules; reaction is then assumed to occur if three of 
the four contact points are correctly matched and within 
a specified maximum distance. The rate constants are 
again found to be about 2 orders of magnitude higher 
than expected from a simple geometric argument, but as 
the approach is not analytical, the result is not readily 
generalizablc. 

In the following, we present a general expression for 
the partially diffusion-controlled rate constant fcnc for 
two spherical molecules with fully asymmetric binding 
patches. The theoretical derivation is given in Sec. ITT1 
Various aspects of our general expression are investi- 
gated in Sec. IIIII where we demonstrate that previous 
approaches are, as expected, limiting cases of our general 
treatment (Sec. lIII All , discuss the dependence of the rate 
constant on orientational constraints fSection lHIBfl . and 
compare numerical values obtained from our expression 
with the result of a Brownian Dynamics simulation by 
Northrup and Erickson CSec. IIII C|) . 



II. THEORY 
A. Model and coordinate system 

Our model for bimolecular association (see Fig.^) con- 
sists of two spherical molecules with radii R\ and i? 2 , 
respectively, whose relative distance and angular orien- 
tation change by translational and rotational diffusion 
with diffusion constants D = D\ lans + D| rans , D\ ot and 
DJ, * . The center of sphere 1 coincides with the origin of 
a fixed-space coordinate system {x, y, z}. The position of 
the center of sphere 2 is specified by the center-to-center 
vector r whose spherical coordinates with respect to the 
fixed-space coordinate system are given by (r, 9, <j>). 

Each sphere carries a body-fixed coordinate system, 
denoted by {x\, yi, z\ \ and {x2, J/2, ^2}, respectively, with 
the axes z\ and Z2 pointing along r when the two spheres 
are perfectly aligned (and hence z\ and z% can be thought 
of pointing at the "center" of the reactive patch). The 
orientation of these body-fixed coordinate systems with 
respect to the fixed-space coordinate system {x,y,z} is 
parametrized by sets of Euler angles 'J'i = (<^i,0i,Xi) 
and ^2 = (4>2, 62, X?)- The angles 4>i and Oi, i = 1,2, are 
the usual azimuthal and polar coordinates of the Zi axis, 
whereas Xi measures the angle from the line of nodes, de- 
fined to be the intersection of the xy and the Xiyi planes, 
to the yt axis. The set (r, 9, 4>, ^1, ^2) comprises the ab- 
solute coordinates of the system. 

For a convenient formulation of the reaction condi- 
tion, we additionally introduce a relative coordinate sys- 




FIG. 1: Absolute and relative coordinate system describing 
the diffusional motion of the two spheres (see text). For the 
sake of clarity, all \ and some of the (p angles have been omit- 
ted in the drawing. 

tern {x rc i, y xe \, z lc \}. The z rc \ axis coincides with the 
center-to-center vector r, whereas the lies in the 

plane spanned by r and the z axis of the fixed-space 
coordinate system {x, y,z}. The Euler angles ^a = 
((/)a,0a,Xa) and ^ B = (</>b,#b,Xb) specify the orien- 
tation of the body-fixed coordinate systems {x\,y\, z\} 
and {x2,y2, Z2] with respect to the coordinate system 

B. Reaction condition 

To fully specify the position and orientation of two 
rigid bodies, nine variables are required, for instance, 
as introduced through our absolute coordinate system, 
(r, 6, (j>, ^ijNf'a). However, for the expression of our re- 
action condition, only five variables, describing the dis- 
tance between the two spheres and their relative orienta- 
tion, are needed (see Fig. |2J). First, the center-to-center 
distance is parametrized by r. The differences in the 
orientation of the two spheres can be fully captured by 
the differences in the Euler angles ^a = {4>a,0a,Xa) 
and <3> B = (4>b,0 b ,Xb), namely, 59 = \9 A - 9 B \, 6<j> = 
\<Pa - 4>b\, and S\ = \xa - Xb\- Finally, we need a 
measure for the extent to which the reactive patches on 
the spheres are aligned with the center-to-center vector 
r, which can be represented by the sum of the polar an- 
gles 9a + Ob- To facilitate the subsequent calculations, 
we replace the conditions on 9a — 9b and 9a + 9b with 
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Sphere 2 




where 



Sphere 1 



FIG. 2: The axes and angles relevant to the reaction condi- 
tion, Eqs. (0. The angles 9a and 9b measure how close the 
center of each reactive patch (coninciding with the body-fixed 
axes zi and Z2, respectively) is to the center-to-center vector 
r. The angles 5<j> and Sx denote relative torsion angles of the 
two body- fixed coordinate systems (a;i, yi, zi) and (x2, 2/2,-22). 
To facilitate the visualization of these two angles, the origin 
of the x\ and yi axes (belonging to the coordinate system of 
sphere 1) has been shifted such as to coincide with the origin 
of the coordinate system of sphere 2. Our reaction condition, 
Eqs. 0, requires near-optimal alignment, i.e., all angles 6a, 
9b, 8(j>, and 5x must be below given limits. 
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is the Laplace operator acting on the center-to-center vec- 
tor r, expressed in the spherical coordinates (r, 9, <fi), and 
the 8 Si , s = x,y,z and i — 1,2, denote an infinitesimal 
rotation of sphere i about its body-fixed Si axis. Equa- 
tion (J2J) can be viewed as composed of three individual 
diffusional contributions, namely, the diffusional motion 
of the center of mass of sphere 2 relative to sphere 1 and 
the rotational diffusion of each sphere. 

As in quantum mechanics, we can define angular mo- 
mentum operators J Si = —ihd/8S Si as generators of in- 
finitesimal rotations of the spheres about their body-fixed 
axes, and can hence rewrite Eq. J2J) as 

= DV 2 c + D[ ot jfc + D r 2 ot J 2 c, (4) 



where Jf = (-ih) 2 {Jf i 
relations 



J 2 

J Vi 



J 2 .). Using the basic 



dS Xi = d9i sin x% ~ <H>i sin 9i cos \i 
dS yi = d9i cos Xi + d<f)i sin 9i sin Xi 
d5 Zi = d(f>i cos 9i + dxi, 



(5) 



we can express the operator J 2 in terms of the Euler 
angles = (cfo, 0i= Xi), 



independent constraints on 9a and 9b- Our reaction con- 
dition is therefore 



r = R 1 +R 2 = R, 
8a.b < 8\bi 



< 



(1) 



M = \xa - Xb\ < Sxo- 



C. Derivation of the rate constant expression 

To derive an expression for the association rate con- 
stant, we determine the concentration c(r, 9, 4>, ^i, ^2) 
of spheres 2 as solution of the steady-state translational- 
rotational diffusion equation in the absolute coordi- 
nate system {r, 9, <fi, \J>i, VE^} introduced in the preceding 
Sec. EH 



dc 
dt 
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The advantage of the formulation of the diffusion equa- 
tion, Eq. ©, in terms of the operators Jf in Eq. |0J, 
lies in the fact that the properties of the Jf are well- 
known, in particular their eigenfunctions, which are 
given by the Wigner rotation matrices T> l mn ((f), 9, x) = 
e- im ^d l mn (9)e- m ^ [|. 

The general solution to Eq. (0} that obeys the bound- 
ary condition at r — > 00, 



lim c(r, 9, <f>, ^1,^2) = Co = const. 



(7) 



can therefore be written as a series of products of the 
eigenfunctions of V 2 -, Jf, and Jf , 



c(r,0,0,*i,tf 2 ) = co + J2 



jmmi)iim2ti2 
/ , / , A llxh 
ll±l 2 mmim 2 nin 2 



where 



flhh(r) 



Ki + i, 2 (£r) 
(fr)V2 



(9) 
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are the modified Bessel functions of the third kind [T(i| 
(giving the desired behavior fu t i 2 (r) — > as r — * oo), 
with £ = [{D\ ot /D)h(h + 1) + (D 2 ot /D)l 2 (l 2 + 1)]V2. 

For the boundary condition at r = _R, the usual, but 
analytically hardly tractable radiation boundary condi- 
tion is 

9c K 
dr „ = D 



F(* A) *B)c(i2,0,^*i,*2), (10) 



where re quantifies the extent of diffusion control in the re- 
action, and A , y B )=H(6° A - 6 A )H(9° B - 9 B )H(6<j) - 
S(f>)H(Sxo~^x) represents the reaction condition Eq. @J, 
where H(x) is the step function defined by TC(x) = for 
x < and H(x) = 1 for x > 0. 

In our approach, we express the radiation boundary 
condition using the constant-flux approximation as in- 
troduced by Shoup et al. |5j , by requiring that the flux is 
a constant over the angular ranges in which the reaction 
can take place, 



dc 
dr 



QF(^ A ,^ B ), 



(11) 



and that Eq. 1)11) ll is obeyed on the average over the sur- 
faces of the spheres, that is, 



= ^JdnJd-i> 1 J dv 2 F(# A , #u) c(R, e, 0, * 2 ), 

(12) 

where we have introduced the abbreviation J <if2 = 
/ sin OdO J d<j>. 

To proceed, we express F{^a, ^b) in absolute coordi- 
nates. First, we expand F(^a, *S?b) in terms of rotation 
matrices, 



F^a^b) 
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V 1 * A*a)V 1 °{* b ), (13) 



where the expansion coefficients (j m ^ nAmBni ' are given 
by 



^mAnAmsne _ + 1 2Zg + 1 



8tt 2 



8tt 2 



rf*; 



x 2^1 a)V 1 °* (* b )F(* a ,*b) 



m A n A v 

+ 1 2/, 



1 
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8tt 2 



(14) 



The absolute coordinate system z} can be trans- 

formed into the relative coordinate system {x Te i, y Ie h z roi} 
by rotations through the three Euler angles (0 — ir, 8, 0). 
The corresponding transformations of the rotation ma- 
trices appearing in Eq. I|13|) are then 

^ A n A = E V ^m A (0 " 0, 0)P^ inA (fx), 

jni 

^ BnB = E V ^ mB fa ~ ^> 0)2^ 2 „ B (M> 2 ). 



The expansion coefficients A 



mminim2n2 
lllh 



in Eq. © 



can be obtained by substituting the expansion for 
F( 1 S r J 4, ^s), Eq. I|13[l . expressed in absolute coordinates 
(r, #,</>, ^1, Vl^), using the above transformations, into 
Eq. ffTTf) - which yields 



\mm\Tl\ 7712^2 



Q 



(-1) 



m+mi+m.2 — ni — ni 



x V ^(2TTT)( 1 



1 h i 2 \ 

m —mi —m°2 J 

^ mA - ni (o^-™J> (15) 

m A 



where 



; re ( ^ ^ ^ ) is the Wigner 3-j symbol. Evaluating 
Eq. i|12|) using the expansion coefficients, Eq. I|15|) . yields 



for the constant Q 
D 

c a o x — a 



x 



2Zi + 1 2l 2 



ihi 
+1 



n h i 2 {R) 



8tt 2 8tt 2 



1 -Hi -Hi 
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(16) 



n— — 2i ra— — l\ 

where we have introduced 

a = J dil J d&t J d* 2J F(*A,*B) 

= (4^) 3 ^ <5 Xo (l-cos^)(l-cos^). (17) 

ao I (47r x 8-7T 2 x 87r 2 ) is the fraction of angular orientational 
space over which the reaction can occur. In derivin 
Eqs. lfT5|> and (fTH|l . we have made use of the identities 



2Z + 1 

47T 



8tt 2 
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The diffusion-controlled rate constant is given by 



feDC 



1 R 2 D 
Bvr 2 ) 2 c 
1 i? 2 L> 

B7T 2 ) 2 " Co" 



,H1 j a'*, / d* 2 ^ 



-a Q. 



(18) 
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Since the functions fu^ir), denned in Eq. JJjJ, obey the 
recursion relation 

flhi 2 (r) = -fihh{r) - U{i+i)hi 2 {r), (19) 

the final expression for the diffusion-limited rate con- 
stant, Eq. (|18[l. becomes 



«DC 



D(Ra /8n 2 ) 2 x 



D 

—a 

K 

K l+i/2 (C 



iKi +1/2 {e)-tL*K l+V2 {t*) 



,2/i + 12^2 + l 



Ef p, mn (I h l 2 \ 



n— — l\ x m——l 



i -1 



(20) 



with C = £R. 



III. RESULTS 



A. Limiting cases 



1. Axially symmetric reactive patches 



S\o = 7r and 9 B = it), we arrive at the model introduced 



by Sole and Stockmayer 0. Then, since 9 B — n and 
J* sin9d9Pi(cos9) = if I ^ 0, only the term l 2 = 
(and hence I = l\) gives a nonzero contribution to the 
sum in Eq. J2IJ. Then Eq. (J2H> reduces to 



fc DC = 2nDR 2 (l -cose° A ) 2 x 



-(1-cos^) 



(/ + i/2)^ +1/2 (e) 
^+i/2(e)-e^+3/2(e) 



sin OaCLOa Ph (cos 0a) 



(22) 



where now £* = R[(D\ ot / D)l(l + 1)} 1/2 , which coincides 
with the result of Sole and Stockmayer Q and Shoup 
e£ a/. @. 

Assuming both spheres to be uniformly reactive, 9 A = 
9 B = 54>q = Sxo = only the term I = l± = 
^2 = contributes, and hence £* = 0. Because 
K 1/2 (n/CK 3/2 (C) -» 1 as T - 0, Eq. (23 be- 
comes, in the fully diffusion-controlled case (k — > oo), 
fee = 4:irDR, which is just the classic Smoluchowski 
diffusion-controlled rate constant for two uniformly reac- 
tive spheres. 



Zhou |6j presented an analytical expression for the as- 
sociation rate constant of two spherical molecules bearing 
axially symmetric patches. In the notation of our model, 
this corresponds to setting 5(f>o = S\o = ^, which makes 
C^H = in Eqs. (O and JSUJ), unless m = n = 0. Using 

PqoOM.x) = 4a(4Ax) = Pi(cos9), where Pi(cos9) 
are the Legendre polynomials, the expression for the rate 
constant, Eq. I|20|) . becomes 

fc DC = 47rL>i? 2 (l-cos6» A ) 2 (I-cos^) 2 

4— (1 -cos0 A )(l -cos6>^) 



-R 



\ " 



^ 2 ^ + i/ 2 (e)-e^ + 3/ 2 (e) 



x (2i + l)(2Ji + l)(2I 3 + l) 



sin Qa^a Ph (cos 9a) 



sm9 B d9 B P lB (cos9 B ) l § ) 



which agrees with the solution presented by Zhou 



(21) 



2. Uniform reactivity 

If we assume that one sphere is uniformly reactive and 
the other has an axially symmetric patch (that is, <5(/>o = 



B. Numerical evaluation 

In the following, we shall assume the reaction to be 
fully diffusion-controlled (k — > oo), and take the radii of 
the two spheres to be identical, R\ — R 2 . Instead of plot- 
ting the absolute value of the association rate constant 
fcoc, we introduce the dimensionless relative association 
rate constant fcjj c = k-Dc/^DR, which is the ratio of 
the orientation-constrained rate constant to the Smolu- 
chowski rate constant for two uniformly reactive spheres. 

The full dependence of the relative association rate 
constant on 9 A , 9 B , S4>o, and <5xo is not easy to display 
in a single plot. For simplicity, we set all four param- 
eters equal, and in Fig. [31 plot the relative association 
rate constant k^ c computed from Eq. (|2Tj|) as a func- 
tion of a single parameter (referred to $o in the follow- 
ing). For comparison, we also show the relative asso- 
ciation rate expected from a purely probabilistic argu- 
ment (geometric rate), given by the fraction of angular 
orientational space over which the reaction can occur, 
a /(47r x 8tt 2 x 8tt 2 ) = $ (1 - cos $o) 2 /4tt 2 . 

It is evident from Fig.[3]that the difference between the 
rate constant fcjj c and the geometric rate gets more strik- 
ing as the angular constraint $o becomes more stringent. 
For instance, in the important case of sterically highly 
specific protein-protein interactions where $o will typ- 
ically range between 5° and 15°, the geometric rate is 
about 2 to 3 orders of magnitude too low, as compared 
with the association rate computed from Eq. (|20|l . 
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FIG. 3: Diffusion-controlled (n — ► 00) relative associa- 
tion rate constant k^c = koc/^DR, with koc computed 
from Eq. I|20|l . as a function of the angular constraint $0 = 
@a — @b — 8<t>o = ^Xo (solid curve). Also shown is the 
rate expected from a simple probabilistic argument, kuc = 
$o(l — cos "I>o) 2 /47r 2 (geometric rate; dashed curve). 



C. Comparison against Brownian dynamics 
simulations 

In the Brownian dynamics simulations by Northrup 
and Erickson protein molecules are modeled as hard 
spheres of R = 18 A diffusing in water (77 ~ 8.9 x 
10~ 4 Ns/m 2 ) at T — 298 K; no forces are assumed to 
act between the molecules. The translational and rota- 
tional diffusion constants are computed from the Stokes- 
Einstein relations L» trans = k B T/6m]R and D Iot = 
kBT/8TTi]R 3 , respectively. 

Instead of angular constraints, the model uses a 
contact-based reaction condition. A set of four distinctly 
numbered contact points is mounted on each sphere in 
a 17 A x 17 A square arrangement on a plane tangen- 
tial to the surface of the sphere. Reaction is assumed 
to occur when at least three of the four contact points 
are correctly matched and within a maximum distance of 
2 A. 

We performed numerical simulations to estimate the 



angles 9 A , 9 B , 6<fro, and 5\o ( as defined in our model, 
see Sec. [HJ that correspond to this contact-based reac- 
tion condition. Clearly, there will be a multiplicity of 
sets of these angles for which the contact-based reaction 
criterion is met. To reduce the search space in a reason- 
able way, we looked for geometric configurations where 
all four angles were equal, 9 A = 9% = 8<p = 5\o, and 
found that the contact-based reaction condition can be 
well represented by an angular constraint of 9 A = 9 B = 

S X 6-7°. 

With these angular constraints, numerical evaluation 
of Eq. (|20|) with the parameters specified above and 
k — > 00 gives /cdc — 1-04 x 10 5 M^ 1 s^ 1 , which is in very 
good agreement with the value obtained from the Brow- 
nian dynamics simulation by Northrup and Erickson |tj , 

foe = l x lf^M" 1 s _1 . 

IV. SUMMARY 

We have presented a general expression for the 
diffusion-controlled association rate of two molecules 
where reaction can occur solely if specified constraints on 
the mutual orientation arc fulfilled. Our solution goes far 
beyond previous treatments in the ability to impose very 
general, asymmetric orientational constraints, as needed 
for instance in a proper description of the sterically highly 
specific association of two proteins. 

Since our expression for the rate constant, Eq. I|2(J|I . 
was derived under the assumption of no forces acting 
between the two molecules, a comparison of measured 
association rates with their theoretical values calculated 
from Eq. H20|l should reveal the extent to which long- 
range interactions contribute to the rate of intermolecular 
association. Such an investigation would be of particular 
interest in the case of the association of proteins with 
small ligands and other proteins. 
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